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ABSTRACT 

The dynamical equations describing the evolution of a self-gravitating fluid can be 
rewritten in the form of a Schrodinger equation coupled to a Poisson equation de- 
termining the gravitational potential. This wave-mechanical representation allows an 
approach to cosmological gravitational instability that has numerous advantages over 
standard fluid-based methods. We explore the usefulness of the Schrodinger approach 
by applying it to a number of simple examples of self-gravitating systems in the weakly 
non-linear regime. We show that consistent description of a cold self-gravitating fluid 
requires an extra "quantum pressure" term to be added to the usual Schrodinger 
equation and we give examples of the effect of this term on the development of grav- 
itational instability. We also show how the simple wave equation can be modified by 
the addition of a non-linear term to incorporate the effects of gas pressure described 
by a polytropic equation-of-state. 
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> . 1 INTRODUCTION 



The local Universe displays a rich hierarchical pattern of 
galaxy clustering that encompasses a large range of length 
scales, culminating in rich clusters and superclusters; the 
early Universe, however, was almost smooth, with only slight 
ripples seen in the cosmic microwave background radia- 
tion (e.g. Smoot et al. 1992). Models of the evolution of 
structure link these observations by appealing to the effect 
of gravitational instability, e.g. Lahav et al. (2002). Low- 
amplitude primordial perturbations of wave-like character 
within a largely homogeneous universe become amplified by 
their own self-gravity, first linearly which then, as their am- 
plitudes build up, results in non-linear collapse. The linear 
theory of perturbation growth for cosmological density fluc- 
tuations is well-established, and is founded on tried-and- 
tested representation of matter on cosmological scales as 
a cold fluid (e.g. Peebles 1980). The non-linear regime is 
much more complicated and generally not amenable to ana- 
lytic solution. Instead, numerical iV-body simulations using 
particles rather than waves to trace cosmic structure, have 
led the way towards an understanding of strongly developed 
clustering. Analytic approaches based on particles, chiefly 
the Zel'dovich (1970) approximation, have also proved use- 
ful. Although these agree to an extent with fully numerical 
calculations (Coles, Melott & Shandarin 1993), they break 
down when particle trajectories cross to form singularities 
known as caustics. 



In this paper we adopt an approach to cosmic structure 
that has advantages over both fluid and particle descrip- 
tions. Following Widrow & Kaiser (1993) and Coles (2002), 
we construct a formalism in which the dynamics of gravita- 
tional is couched in the language of quantum mechanics and 
governed by a Schrodinger wave equation. The advantages 
of this approach will be made apparent as we explore its be- 
haviour in simple examples which can be solved exactly and, 
within the Schodinger approach, to various levels of approx- 
imation. Our aim is to construct a useful tool for evolving 
density fluctuations into the mildly non-linear regime, be- 
yond the breakdown of linear perturbation theory. 

The importance of the search for analytic understand- 
ing of structure formation is often overlooked in the light 
of the tremendous advances that have been made recently 
in computational astrophysics. One motivation for a wider 
range of theoretical techniques is practical. High-resolution 
A-body experiments are CPU-intensive and it is difficult to 
find sufficient resources to use them to explore a large range 
of choices of initial fluctuations, different cosmological pa- 
rameters, and so on. A neat analytic approximation can do 
this job more efficiently, at least to the point of narrow- 
ing down the field of contenders for fuller study. The other 
motivation is methodological: the complex "cosmic web" of 
large-scale structure requires explanation rather than repro- 
duction, and explanation rests on using analytic approaches 
to identify the important physics as much as possible. 
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The layout of this paper is as follows. In Section 2 we 
run through the standard basics of gravitational instabil- 
ity in an expanding universe, including perturbation the- 
ory and the Zel'dovich approximation. We then, in Section 
3, present the details of the Schrodinger approach and de- 
rive a spherically-symmetric solution. In Section 4 we apply 
the method to examples of one-dimensional collapse show- 
ing how various approximations within the wave-mechanical 
framework behave when compared to the real data. We sum 
up in Section 5. 



2 COSMOLOGICAL STRUCTURE 
FORMATION 

2.1 The Background Cosmology 

Space-times consistent with the Cosmological Principle can 
be described by the Robertson-Walker metric 



dr 2 



1 — K7" 



+ r 2 d6 2 + r 2 sm 2 6d(j> 2 ] , (1) 



where k is the spatial curvature, scaled so as to take the 
values or ±1. The case K = represents flat space sec- 
tions, and the other two cases are space sections of con- 
stant positive or negative curvature, respectively. The time 
coordinate t is cosmological proper time and a(t) is the cos- 
mic scale factor. The dynamics of a Friedmann-Robertson- 
Walker universe are determined by the Einstein gravita- 
tional field equations which can be written 



8ttGp ^- + Ac , 
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(2) 
(3) 
(4) 



These equations determine the time evolution of the cosmic 
scale factor a(t) (the dots denote derivatives with respect to 
cosmological proper time i) and therefore describe the global 
expansion or contraction of the universe. The behaviour of 
these models can further be parametrized in terms of the 
Hubble parameter H = (a/ a) and the density parameter 
f2 = 8irGp/3H 2 , a suffix representing the value of these 
quantities at the present epoch when t — to. 



2.2 Fluid Treatment 

In the standard treatment of the Jeans Instability one begins 
with the dynamical equations governing the behaviour of a 
self-gravitating fluid. These are the Euler equation 

d(v) 1 

-^ + (v- V)v + -Vp + V0 = O ; (5) 

the continuity equation 

g + V(pv) = , (6) 

expressing the conservation of matter; and the Poisson equa- 
tion 



describing Newtonian gravity. If the length scale of the per- 
turbations is smaller than the effective cosmological horizon 
dn = c/H, a Newtonian treatment of cosmic structure for- 
mation is still expected to be valid in expanding world mod- 
els. In an expanding cosmological background, the Newto- 
nian equations governing the motion of gravitating particles 
can be written in terms of 



r/a(t) 



(8) 



(the comoving spatial coordinate, which is fixed for observers 
moving with the Hubble expansion), 



v = r — Hr = ax 



0) 



(the peculiar velocity field, representing departures of the 
matter motion from pure Hubble expansion), p(n,t) (the 
matter density), and <f)(x,t) (the peculiar Newtonian gravi- 
tational potential, i.e. the fluctuations in potential with re- 
spect to the homogeneous background) determined by the 
Poisson equation in the form 



V x 4> = 47rGa (p — po) = AnGa po5. 



(10) 



In this equation and the following the suffix on Va, indicates 
derivatives with respect to the new comoving coordinates. 
Here po is the mean background density, and 



po 



(11) 



is the density contrast. Using these variables the Euler equa- 
tion becomes 



-^rr 1 + V ■ V x V = VxP - V x 

Ot p 



(12) 



The first term on the right-hand-side of equation 112H arises 
from pressure gradients, and is neglected in dust-dominated 
cosmologies. Pressure effects may nevertheless be important 
in the the (collisional) baryonic component of the mass dis- 
tribution when nonlinear structures eventually form. The 
second term on the right-hand side of equation 1121 is the 
peculiar gravitational force, which can be written in terms 
of g = — Vx</>/a, the peculiar gravitational acceleration of 
the fluid clement. If the velocity flow is irrotational, v can 
be rewritten in terms of a velocity potential (p v : 



v = -V x <f> v /a 



(13) 



: 4irGp 



(7) 



This is expected to be the case in the cosmological setting 
because (a) there are no sources of vorticity in these equa- 
tions and (b) vortical perturbation modes decay with the 
expansion. We also have a revised form of the continuity 
equation: 

g+3i/p+iVx(pv) = . (14) 



2.3 Linear Perturbation Theory 

The procedure for handling perturbations is to linearise the 
Euler, continuity and Poisson equations by perturbing phys- 
ical quantities defined as functions of Eulerian coordinates, 
i.e. relative to an unperturbed coordinate system. Expand- 
ing p, v and 4> perturbatively and keeping only the first- 
order terms in equations 1121 and 114H gives the linearised 
continuity equation: 
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dt a 



(15) 



which can be inverted, with a suitable choice of boundary 
conditions, to yield 



8 = 



1 



aHf 



(V* 



(16) 



The function / ~ f2o' 6 ; this is simply a fitting formula to 
the full solution (Peebles 1980). The linearised Euler and 
Poisson equations are 

9v 



a 

ot a 



pa 



-V x? 



V x 2 </> = 47rGa 2 poS; 



;V X 0, (17) 

(18) 

|vj, |0|, \S\ < 1 in equations 0, El & GU- From tnese 
equations, and if one ignores pressure forces, it is easy to 
obtain an equation for the evolution of 5: 



8 + 2HS- ^H 2 8 = 0. 



(19) 



For a spatially flat universe dominated by pressureless mat- 
ter, po{t) = l/67rGi 2 and equation 11911 admits two linearly 
independent power law solutions 5(x, t) = D±(t)S(x), where 
D + {t) oc a(t) oc t 2 / 3 is the growing mode and -D-(t) oc t~ 1 
is the decaying mode. 

The above considerations apply to the evolution of a 
single Fourier mode of the density field <5(x,i) = D+(t)<5(x). 
What is more likely to be relevant, however, is the case of a 
superposition of waves, resulting from some kind of stochas- 
tic process in which the density field consists of a superpo- 
sition of such modes with different amplitudes. A statistical 
description of the initial perturbations is therefore required, 
and any comparison between theory and observations will 
also have to be statistical. Many versions of the inflation- 
ary scenario for the very early universe (Guth 1981; Guth 
& Pi 1982; Linde 1982; Albrecht & Steinhardt 1982) predict 
the initial density fluctuations to take the form of a Gaus- 
sian random field in which the initial Fourier modes of the 
perturbation field have random phases. 

The linearised equations of motion provide an excellent 
description of gravitational instability at very early times 
when density fluctuations are still small (S <C 1). The lin- 
ear regime of gravitational instability breaks down when 8 
becomes comparable to unity, marking the commencement 
of the quasi-linear (or weakly non-linear) regime. During 
this regime the density contrast may remain small (8 < 1), 
but the phases of the Fourier components <5k become sub- 
stantially different from their initial values resulting in the 
gradual development of a non-Gaussian distribution func- 
tion even if the primordial density field is Gaussian owing 
to mode-mode interactions. 

Perturbation theory fails when non-linearity develops, 
but the fluid treatment is intrinsically approximate anyway. 
A fuller treatment of the problem requires a solution of the 
Boltzmann equation for the full phase-space distribution of 
the system /(x, v, t) coupled to the Poisson equation Q 
that determines the gravitational potential. In cases where 
the matter component is collisionless, the Boltzmann equa- 
tion takes the form of a Vlasov equation: 



dt ^-^ \ dxi dvi ' dxi J 



(20) 



The fluid approach outline above can only describe cold ma- 
terial where the velocity dispersion of particles is negligible. 
But even if the dark matter is cold, there may be hot com- 
ponents of baryonic material whose behaviour needs also to 
be understood. Moreover, the fluid approach assumes the 
existence of a single fluid velocity at every spatial position. 
It therefore fails when orbits cross and multi-streaming gen- 
erates a range of particle velocities through a given point. 

2.4 The Zel'dovich Approximation 

Fortunately the formation of the main elements of large- 
scale structure can nevertheless be understood using rela- 
tively simple tools, such as the Zel'dovich approximation 
(Zel'dovich 1970). Let the initial (i.e. Lagrangian) coordi- 
nate of a particle in this unperturbed distribution be q. Now 
each particle is subjected to a displacement corresponding 
to a density perturbation. In the Zel'dovich approximation 
the Eulerian coordinate of the particle at time t is 



r(i,q) =a(t)[q-b(t)V q $„(q)], 



(21) 



where r = a(i)x, with x a comoving coordinate, and we have 
made a(t) dimensionless by dividing throughout by a(ti), 
where ti is some reference time which we take to be the 
initial time. Notice that particles in the Zel'dovich approx- 
imation execute a kind of inertial motion on straight line 
trajectories. The derivative on the right hand side is taken 
with respect to the Lagrangian coordinates. The dimension- 
less function b(t) describes the evolution of a perturbation in 
the linear regime, with the condition b(ti) = 0, and therefore 
solves the equation 



6 + 2 (^) b - AnGpb 







(22) 



cf. equation (19). For a flat matter-dominated universe we 
have b oc t 2 ^ 3 as before. The quantity <l , o(q) is proportional 
to a velocity potential, of the type introduced above, i.e. a 
quantity of which the velocity field is the gradient: 

V=--tfr = a— = -a&V,$o(q); (23) 

this means that the velocity field is irrotational. The quan- 
tity $o (q) is related to the density perturbation in the linear 
regime by the relation 



8 = 6Vi*o, 



(24) 



which is a simple consequence of Poisson's equation. 

The Zel'dovich approximation, though simple, has a 
number of interesting properties. First, it is exact for the 
case of one dimensional perturbations up to the moment 
of shell crossing. As we have mentioned above, it also in- 
corporates irrotational motion, which is required to be the 
case if it is generated only by the action of gravity (due 
to the Kelvin circulation theorem) . For small displacements 
between r and a(t)q, one recovers the usual (Eulerian) lin- 
ear regime: in fact, equation (22) defines a unique mapping 
between the coordinates q and r (as long as trajectories do 
not cross); this means that p(r,t)d 3 r — {p(ti))d 3 q or 



PM) = 



\J(r,t)\ 



(25) 



where | J(r , t) | is the determinant of the Jacobian of the map- 
ping between q and r: 9r/9q. Since the flow is irrotational 
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the matrix J is symmetric and can therefore be locally di- 
agonalised. Hence 

3 

p(r,t) = (p(t)>JJ[l + 6(t)o i (q)]- 1 : (26) 

£=1 

the quantities 1 + b(t)oti are the eigenvalues of the matrix J 
(the ai are the eigenvalues of the deformation tensor). For 
times close to ti, when \b(t) at\ <C 1, equation (26) yields 

5~ -(qi +a 2 +a 3 )b(t), (27) 

which is just the law of perturbation growth obtained by 
solving equation 119H . 

At some time t sc , when b(t 3C ) — —1/ctj, a singularity 
appears and the density becomes formally infinite in a re- 
gion where at least one of the eigenvalues (in this case a.j) is 
negative. This condition corresponds to the situation where 
two points with different Lagrangian coordinates end up at 
the same Eulerian coordinate. In other words, particle tra- 
jectories have crossed and the mapping between Lagrangian 
and Eulerian space is no longer unique. A region where the 
shell-crossing occurs is called a caustic. For a fluid element 
to be collapsing, at least one of the «j must be negative. If 
more than one is negative, then collapse will occur first along 
the axis corresponding to the most negative eigenvalue. If 
there is no special symmetry, one therefore expects collapse 
to be generically one-dimensional, from three dimensions to 
two. We shall make use of the fact that in one-dimensional 
cases, the Zel'dovich solution is exact until the moment of 
shell- crossing. 

The Zel'dovich approximation matches very well the 
evolution of density perturbations in full iV-body calcu- 
lations until the point where shell crossing occurs (Coles, 
Melott & Shandarin 1993). After this, the approximation 
breaks down completely. Particles continue to move through 
the caustic in the same direction as they did before. Parti- 
cles entering a pancake from either side merely sail through 
it and pass out the opposite side. The pancake therefore ap- 
pears only instantaneously and is rapidly smeared out. In 
reality, the matter in the caustic would feel the strong grav- 
ity there and be pulled back towards it before it could escape 
through the other side. Since the Zel'dovich approximation 
is only kinematic it does not account for these close-range 
forces and the behaviour in the strongly non-linear regime is 
therefore described very poorly. Furthermore, this approx- 
imation cannot describe the formation of shocks and phe- 
nomena associated with pressure. Attempts to understand 
properties of non-linear structure using a fluid description 
therefore resort to further approximations (Sahni & Coles 
1995) to extend the approach beyond shell-crossing. One 
method involves the so-called adhesion model (Gurbatov, 
Saichev & Shandarin 1989). In this model one assumes that 
the particles stick to each other when they enter a caustic 
region because of an artificial viscosity which is intended to 
simulate the action of strong gravitational effects inside the 
overdensity forming there. This 'sticking' results in a cancel- 
lation of the component of the velocity of the particle per- 
pendicular to the caustic. If the caustic is two-dimensional, 
the particles will move in its plane until they reach a one- 
dimensional interface between two such planes. This would 
then form a filament. Motion perpendicular to the filament 
would be cancelled, and the particles will flow along it until 



a point where two or more filaments intersect, thus forming 
a node. The smaller is the viscosity term, the thinner will 
be the sheets and filaments, and the more point-like will 
be the nodes. Outside these structures, the Zel'dovich ap- 
proximation is still valid to high accuracy. See Shandarin & 
Zel'dovich (1989) for further discussion. 



2.5 The Strongly Nonlinear Regime 

Further into the non-linear regime, bound structures form. 
The baryonic content of these objects may then become im- 
portant dynamically: hydrodynamical effects (e.g. shocks), 
star formation and heating and cooling of gas all come into 
play. The spatial distribution of galaxies may therefore be 
very different from the distribution of the (dark) matter, 
even on large scales. Attempts are only just being made to 
model some of these processes with cosmological hydrody- 
namics codes, such as those based on Smoothed Particle Hy- 
drodynamics (SPH; Monaghan 1992), but it is some measure 
of the difficulty of understanding the formation of galaxies 
and clusters that most studies have only just begun to at- 
tempt to include modelling the detailed physics of galaxy 
formation. In the front rank of theoretical efforts in this 
area are the so-called semi-analytical models which encode 
simple rules for the formation of stars within a framework 
of merger trees that allows the hierarchical nature of gravi- 
tational instability to be explicitly taken into account. 

The spatial distribution of particles obtained using the 
adhesion approximation represents a sort of "skeleton" of 
the real structure: non-linear evolution generically leads to 
the formation of a quasi-cellular structure, which is a kind 
of "tessellation" of irregular polyhedra having pancakes for 
faces, filaments for edges and nodes at the vertices . This 
skeleton, however, evolves continuously as structures merge 
and disrupt each other through tidal forces; gradually, as 
evolution proceeds, the characteristic scale of the structures 
increases. In order to interpret the observations we have al- 
ready described, one can think of the giant "voids" as being 
the regions internal to the cells, while the cell nodes corre- 
spond to giant clusters of galaxies. While analytical meth- 
ods, such as the adhesion model, are useful for mapping 
out the skeleton of structure formed during the non-linear 
phase, they are not adequate for describing the highly non- 
linear evolution within the densest clusters and superclus- 
ters. In particular, the adhesion model cannot be used to 
treat the process of merging and fragmentation of pancakes 
and filaments due to their own (local) gravitational instabil- 
ities, which must be done using full numerical computations. 



3 WAVE MECHANICS AND STRUCTURE 
FORMATION 

3.1 Introduction 

The combination of Eulerian fluid-based perturbation the- 
ory, Lagrangian approaches like the Zel'dovich approxima- 
tion and brute-force iV-body numerics has been successful 
at establishing a standard framework within which the basic 
features of large-scale structure can be described and under- 
stood. There remain significant obstacles to fuller analytic 
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description. The three most interesting from the perspective 
followed in this paper are: 

(i) Perturbation Theory. Standard perturbation meth- 
ods do not guarantee a density field that is everywhere pos- 
itive. This basic problem is obvious in the standard case 
where there is a Gaussian distribution of initial fluctuations. 
When the variance of the distribution of 5 is of order unity 
(its mean is, by definition, always unity) then a Gaussian 
distribution assigns a non-zero probability to regions with 

5 < — 1, i.e. with p < 0. 

(ii) Shell-crossing. Methods such as the Zel'dovich ap- 
proximation and its variations break down at shell-crossing, 
where the density becomes infinite. 

(iii) Gas Pressure. Analytic techniques for modelling 
the effects of gas pressure are scarce, even in relatively simple 
systems such as the Lyman-a absorbing systems (Matarrese 

6 Mohayee 2002). Fully numerical techniques, such as SPH, 
are generally expected to be the way forward but even they 
are not yet able to incorporate all relevant gas physics. 



3.2 The Widrow— Kaiser Approach 

A novel approach to the study of collisionless matter, 
with applications to structure formation, was suggested by 
Widrow & Kaiser (1993). It involves re- writing of the fluid 
equations given in Section 2 in the form of a non-linear 
Schrodinger equation. The equivalence between this and the 
fluid approach has been known for some time; see Spiegel 
(1980) for historical comments. Originally the interest was 
to find a fluid interpretation of quantum mechanical effects, 
but in this context we shall use it to describe an entirely 
classical system. 

To begin with, consider the continuity equation and 
Euler equation for a curl-free flow in which v = V0), in 
response to some general potential V . In this case the con- 
tinuity equation can be written 

g + V • (pV</>) = . (28) 

It is convenient to take the first integral of the Euler equa- 
tion, giving 

9 i^m 2 = -V, (29) 

which is usually known as the Bernoulli equation. The trick 
then is to make a transformation of the form 



V 2 $ = A-KGtpip* 



(33) 



tp — aexp(i<f>/v) 



(30) 



where p = tpip* = \tp\ 2 = a 2 ; the wavefunction tp(x.,t) ev- 
idently complex. Notice that the dimensions of v are the 
same as <j>, namely [L 2 T -1 ]. After some algebra it emerges 
that the two equations above can be written in one equation 
of the form 



where 

„ V 2 V 2 Q 



(31) 



(32) 



To accommodate gravity we need to couple equation 131H to 
the Poisson equation in the form 



taking V to be $. 

This system looks very similar to a Schrodinger equa- 
tion, except for the extra "operator" P. The role of this 
term becomes clearer if one leaves it out of equation 13H 
and works backwards. The result is an extra term in the 
Bernouilli equation that resembles a pressure gradient. This 
is often called the "quantum pressure" that arises when one 
tries to understand a quantum system in terms of a classical 
fluid behaviour. Leaving it out to model a collisionless fluid 
can be justified only if a varies only slowly on the scales 
of interest. On the other hand one can model situations in 
which one wishes to model genuine effects of pressure by ad- 
justing (or omitting) this term in the wave equation. Widrow 
& Kaiser (1993) advocated simply writing 



.. dtp U _ 2 , , f \ i 



(34) 



i.e. simply ignoring the quantum pressure term. In this 
spirit, the constant Ti is taken to be an adjustable parameter 
that controls the spatial resolution A through a de Brogiie 
relation A = h/mv. In terms of the parameter v used above, 
we have v — h/m, giving the correct dimensions for Planck's 
constant. Note that the wavefunction tp encodes the velocity 
part of phase space in its argument through the ansatz 



V'(x) = y/p{x)exp[i9(x)/h], 



(35) 



where V6>(x) = p(x), the local 'momentum field'. This for- 
malism thus yields an elegant description of both the density 
and velocity fields in a single function. 

It should be clear how this approach bypasses imme- 
diately the first two items of the list given in Section 3.1. 
First, the construction 13011 requires that p = \tp\ 2 is positive. 
Secondly, since particles are not treated as point-like enti- 
ties with definite trajectories, shell-crossing does not have 
the catastrophic from in this approach that it does in the 
Zel'dovich approximation. Note that no singularities occur 
in the wavefunction at any time. Although the simplest 
ansatz 1351 does assume a unique velocity at every fluid 
location, it is possible to construct more complex represen- 
tations that allow for multi-streaming (Widrow & Kaiser 
1993) although we shall not use them in this paper. 



3.3 A Spherically-Symmetric Solution 

In order to illustrate how the wave equation approach maps 
onto a fluid description, it is instructive to examine a famil- 
iar problem: a static spherically symmetric self-gravitating 
cloud in hydrostatic balance. This also allows us to explain 
how item (3) of the list in Section 3.1 can be addressed. We 
model the pressure gradients needed to hold an object in 
equilibrium against its self-gravity by the addition to the po- 
tential of a term of the form « 2 1 2 (with k an appropriately- 
chosen constant), which converts the original equation l|34ft 
into a nonlinear Schrodinger equation: 

iv^- = -^V 2 ^ + $(x)ip + n 2 \i(j\ 2 tp + Pi(j (36) 

(Sulem & Sulem 1999), where the final term is the quantum 
pressure. We will now look for solutions of the form 



V , = i? exp(^) 



(37) 
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but we are looking for a static solution so obviously the fluid 
flow velocity should be zero. We therefore need to set (f> to 
be constant in the above transformation so that V0 can be 
zero everywhere. Without loss of generality we can exploit 
the global U(l) symmetry of quantum mechanics to set the 
phase to zero. Hence we can set ip = R and consequently 



\7 2 ip = V 2 -R. The quantum pressure term has the form 



P = 



V 2 R 



R 



(38) 



so that this now cancels the first term on the right-hand- 
side of equation 1361 . In this event, and requiring a static 
solution for ip we get 



$V + K 2 \tp\ 2 1p = 



(39) 



to be solved simultaneous with Poisson's equation. We there- 
fore obtain a simple Helmholtz equation for $: 



(40) 



The solution for spherically symmetry is straightforward: 



$ = $, 



sin Ar 
Ar 



where 

A 2 = — 

K 2 



(41) 



(42) 



and <E> C is constant. There is a particularly neat relation for 
p using equation il 1IH yielding 



4ttGp+-^-<!> = 0, 
which gives 



p(r) = p, 



sin Ar 
Ar ' 



(43) 



(44) 



The resulting density profile will be recognized as that of 
an n = 1 polytrope, i.e. a fluid in which the gas pressure is 
related to the density via p = Kp 2 , with K a constant. 

The important point emerging from this analysis is that 
one can only obtain a static solution in the presence of the 
quantum pressure term. A self-consistent solution requires 
this term to be included and its effect is to stabilize the 
system. In order to apply the wave mechanical approach to 
structure formation, therefore, one has to be more careful 
with this term than Widrow & Kaiser (1993) in their equa- 
tion 13411 . We investigate this issue further in the following 
Section. 



4 ONE-DIMENSIONAL COLLAPSE 

In this section we look at some simple one-dimensional ex- 
amples of gravitational instability using the Schrodinger 
equation. We do this by comparing the solutions we ob- 
tain with the Zel'dovich approximation, which gives an ex- 
act solution up to the point of shell-crossing (Shandarin & 
Zel'dovich 1989). We also want to consider the effects of 
varying the parameter v in the Schrodinger equation and 
how this changes the spatial and velocity resolution of the 
solution, as well as elucidating more clearly the effects the 
"quantum pressure" term. 



4.1 The "Free-Particle" Solution to the 
Schrodinger Equation 

The equation for the trajectory of a particle using the 
Zel'dovich approximation is identical in form to that of iner- 
tial (or ballistic) motion, with the displace replaced by the 
comoving co-ordinate x and the time replaced by b as defined 
in equation (22). The most obviously analogous approach to 
take with the Schrodinger equation is obtained be setting 
the potential V equal to zero, replacing t with b and taking 
spatial derivatives with respect to the comoving co-ordinate 
x. Equation (31) then becomes 



. dib v 2 _ 2 , 



(45) 



or, in one-dimension, 
. dip _ v 2 d 2 i> 

w ~db ~ ~YaP" 

This is the "free-particle" Schrodinger equation and the so- 
lutions are given by 



(46) 



ip(x,t) 



E 

k 



A k e (kx -" b) , 



(47) 



with u> = k 2 v/2. The boundary conditions are given by the 
initial value of the wave function at b = 0: 



ip(x,0) = ^A k e 



(48) 



The values of Ak can be calculated using the Fast Fourier 
Transform of the initial wavefunction, and the value of the 
wavefunction at any future time can be calculated directly 
using equation (47). 

In order to show the equivalence between the Zel'dovich 
approximation and Schrodinger method we consider the two 
solutions when the initial velocity potential is given by a 
single sine wave. In the Schrodinger equation we start with 
an initial wave function given by 



tp(x,0) — - v /poexp(— icosx/v). 



(49) 



This is equivalent to setting the initial density contrast 5 — 
and the velocity v = sin a;. In one dimension, the Zel'dovich 
approximation for the position of a particle at time t is given 
by, 



x{q,t) = q-b 



d$o(g) 



(50) 



cf. equation (21). Conservation of mass is enforced by 
podq = pdx (51) 
and therefore the density contrast is given by, 

,d 2 $ ( q y~' 



i 



i 



(52) 



If the velocity potential 4>o{q) = cosq then the density con- 
trast at time t is then given by 

5 = (1 + bcosx)' 1 - 1 (53) 

and shell-crossing occurs when b — 1 when the density con- 
trast becomes infinite. Figure 1 shows the density contrast 
obtained using both of these methods at three different times 
and for three different values of v. When b = 0.3 density 
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Figure 1. Comparison of the Zel'dovich approximation to the "free-particle" Schrodinger equation. The Zel'dovich approximation is 
given by the dashed line and the solution to the Schrodinger equation by the solid line. At early times the two methods are almost 
identical. Close to shell crossing the results are similar, but with the Schrodinger result behaving more smoothly than the Zel'dovich 
solution. The effect of varying v can be most clearly seen after shell-crossing. 



growth is still approximately linear and the density contrast 
given by the Schrodinger and Zel'dovich methods are almost 
identical. When b — 0.9, just before shell-crossing, the den- 
sity growth has become highly non-linear. The very high 
density peaks seen in the Zel'dovich approximation are not 
seen in the Shroodinger solution. 

The value of v corresponds to a "de Broglie wavelength" 
which gives the spatial resolution of the solution. The ef- 
fects of this can be seen most clearly when b — 2.0 after 
shell-crossing has occurred. Although the Schrodinger and 
Zel'dovich methods show qualitatively different behaviour 
after shell-crossing, this is not physical, and the overall be- 
haviour of the results is in any case similar, particularly for 
small values of v. 

4.2 The Schrodinger equation with 
time-independent potential 

Setting the potential V(x) — in the Schrodinger equation 
gives a good agreement with the Zel'dovich approximation. 
However, in both methods the particles continue on their 
trajectories after shell-crossing with no reference to the self- 
gravity of the collapsing system. Here we look at how the 
inclusion of a potential term in the Schrodinger equation 



might give a more realistic solution after shell-crossing. This 
approach is motivated by the realisation that, at least dur- 
ing linear perturbation growth, the gravitational potential 
corresponding to a density fluctuation barely changes. In a 
different context this has motivated the so-called "frozen- 
potential" approximation (Matarrese et al. 1992; Brainerd, 
Scherrer & Villumsen 1993; Bagla & Padmanabhan 1994), 
which is similar in spirit to the approach we adopt here. 

The simplest way to include the potential is to include 
its initial value, but not to allow this to vary with time. For 
this example we again ignore the quantum pressure term. 
We therefore need to solve the equation 

• 9ip v 2 d 2 i> 

w ¥ = -2^ +cosW (54) 

It should be noted that in this example no account is taken 
of the expansion of the universe. This means that a direct 
comparison with this Zel'dovich approximation is not as 
straightforward as in the previous example and in partic- 
ular the time coordinate is not the same in the two meth- 
ods. However a comparison may be made by considering the 
growth at early times in the linear regime. In the Zel'dovich 
approximation the density contrast is given by 



8 Coles & Spencer 



e— dependent potential 




nculding quantum pressure 







Figure 2. Comparison of the Zel'dovich approximation to the Schrodinger equation with v = 0.1. In the first column the potential V 
does not vary with time. In the second column potential is calculated according to the Poisson equation and in the third column the 
quantum pressure term is included. 



8~b 



dx 2 



-bcosx, 



(55) 



(56) 



and the velocity in comoving coordinates is given by 

dx i^o ; . 
v = — ~ — b— — ~ dsrni. 
at ax 

In a non-expanding universe with a constant potential given 
by V(x) = Vo cosx the density contrast at early times is 

S~-^-t 2 cosx. (57) 

We can therefore compare the Zel'dovich approxima tion at 
time b with the Schrodinger equation at time t = y // 2b/Vo- 
Starting with the initial wavefunction given by 



ip(x,0) = ,/po, 



(58) 



and evolving the wavefunction according to equation (54) 
using Cayley's method (Goldberg, Schcy & Schwartz 1967) 
leads to the results shown in the first column of Figure 2. 
These are for v = 0.1 at times corresponding to the values 
of b = 0.3, 0.9 and 2.0 used before. Comparison with Figure 
1 shows that the extensive smearing of density peaks which 
occurs after shell-crossing in both the "free-particle" and 
Zel'dovich solution is no longer seen in the case where the 
initial potential is included. 



4.3 The Schrodinger equation with time-varying 
potential 

In the second column of Figure 2 we show the solution to 
the Schrodinger equation with the full potential calculated 
using the Poisson equation (33). We start be considering 
the equation without the quantum pressure term. In this 
case we need to solve the coupled Schrodinger and Poisson 
equations: 



. dip v 2 d 2 ip „, , , 



d 2 v 

dx 2 



(59) 



(60) 



To set up the initial wavefunction we choose a single sine 
wave for the density profile, 



S — —S cos x. 



(61) 



In the linear regime the density contrast at time t is given 
by 



S ~ — <$o exp(v / 47rGt) cos x 

and the velocity is 

v — V AuGSq exp(V 4-7rGt) sin x. 



(62) 



(63) 
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Accordingly we compare the Zel'dovich approximation 
at time b with the Schrodinger equation at time t = 
l/V47rGm(&/<5o)- We start with the initial wavefunc- 
tion, using equation (35), given by, 

4>(x,0) — \J po(l — So cos a;) cxp (— WA-kGSq coax/v) , (64) 

and evolve the wavefunction according to equation (59). In 
this case the potential V(x) is calculated by solving the Pois- 
son equation numerically at each time step. The results are 
shown in the second column of Figure 2 for v = 0.1 at times 
corresponding to b — 0.3, 0.9 and 2.0 used before. 

The third column of Figure 2 shows the results for the 
same initial wavefunction, but this time including the quan- 
tum pressure term, as in equation (31). The inclusion of this 
term makes a significant difference to the results at later 
times, particularly after shell-crossing, but is slight earlier 
on during the evolution. This means that useful approxima- 
tions for the mildly non-linear regime can still be be obtained 
even if this term is ignored. 

4.4 Solutions in an expanding background 

Although we do not make further use of it in this paper, 
it is useful to end this section with a few comments about 
what arises when one places the system in an expanding 
background. To see what happens, let us define a scaled 
density x = p/po = (1 + S) and take Q, = 1. The continuity 
equation then becomes 

^ + V • ( X V0) = , (65) 

where the velocity potential 4> is now such that u = dx/dt — 
V(j> and a is the scale factor. It is convenient to take the first 
integral of the Euler equation, giving 

f a+ \ m f = -l^0), (66) 

where 8 = 2&/3a 3 H 2 and <fr is the gravitational poten- 
tial. After some more algebra the system again becomes a 
Schrodingcr-like wave equation, but in a rather than in t 
and using ip 2 = x, such that 

iv llh = -yV 2 V' + , l / V + ^ , (67) 

with V = 4> + 9 and P as before. It is relatively straightfor- 
ward to incorporate the extra terms. 



5 DISCUSSION 

In this paper we have set out an approach to the forma- 
tion of cosmic structure from evolving cosmological density 
fluctuations that relies on a transformation of the evolution 
equations into a Schrodinger- Poisson system. One advantage 
of this new formalism was pointed out by Coles (2002): it 
yields a rather more convincing approach to understanding 
the origin of spatial intermittency in the galaxy distribution 
than that offered by Jones (1999). It also makes a connection 
in the underlying physics with other systems that display 
similar phenomena. In this paper we have argued further 
that it provides an analytical approach that offers benefits 
over standard approaches of both Lagrangian and Eulerian 



types. On the one hand, it always guarantees a positive den- 
sity and on the other it does not produce singularities at 
shell- crossing events. It is clearly important to investigate 
both of these issues further. Szapudi & Kaiser (2003) have 
already completed a study of perturbation theory and its 
statistical ramifications within this approximation. Dealing 
fully with post-shell-crossing phenomena is beyond the scope 
of this paper as it requires a more sophisticated representa- 
tion of the wave-function in velocity space; see Widrow & 
Kaiser (1993) for discussion. 

In this paper we showed how this approach works in 
practice by applying it to a static spherically-symmetric case 
and to some examples of one-dimensional collapse. In doing 
this we looked at a "free-particle" solution of the system that 
resembles the Zel'dovich approximation in some respects. In- 
corporating a fixed gravitational potential yields a solution 
that has some similarities to the "frozen-potential" approx- 
imation (Matarrese et al. 1992). We think the Schrodinger 
formalism provides a useful tool for modelling fluctuations 
in the mildly non-linear regime. In future work we will look 
at its use in modelling the intergalactic medium in order to 
understand properties of Lyman-a absorption systems. 
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